#ifndef GENERATOR_POLYA_HEADER
#define GENERATOR_POLYA_HEADER

#include<string>
#include<fstream>
#include<sstream>
#include<iostream>
#include<algorithm>
#include<iterator>
#include<vector>
#include<map>
#include<cstdlib>
#include<ctime>
#include<cmath>
#include <gsl/gsl_randist.h>
#include "def.hpp"

// gsl_ran_dirichlet (const gsl_rng * r, size_t K, const double alpha[], double theta[])
// gsl_ran_multinomial (const gsl_rng * r, size_t K, unsigned int N, const double p[], unsigned int n[])

using namespace std;

int GE_Polya(gsl_rng* r){

  ofstream outFile;
  ofstream logFile;
  string out_name = std::string(PuzaDEF::Instance()->output_file_name);
  string property_name = std::string(PuzaDEF::Instance()->output_file_name);
  out_name.append(".dat");
  property_name.append(".properties");
  outFile.open(out_name.c_str());
  logFile.open(property_name.c_str());
  
  double* temp_theta = new double[PuzaDEF::Instance()->GE_Polya_K];
  double* temp_alpha = new double[PuzaDEF::Instance()->GE_Polya_K];
  
  for(int k=0; k < PuzaDEF::Instance()->GE_Polya_K; k++){
    temp_alpha[k] = PuzaDEF::Instance()->GE_Polya_Alpha * PuzaDEF::Instance()->GE_Polya_Base[k];
  }


  for(int i =0; i < PuzaDEF::Instance()->GE_Samples; i++){
    gsl_ran_dirichlet(r, PuzaDEF::Instance()->GE_Polya_K, temp_alpha, temp_theta);
    gsl_ran_multinomial(r, PuzaDEF::Instance()->GE_Polya_K, PuzaDEF::Instance()->GE_Polya_N, temp_theta, PuzaDEF::Instance()->GE_Polya_Data);

    outFile << i << "\t";
    for(int k=0; k < PuzaDEF::Instance()->GE_Polya_K; k++){
      outFile << k << ":" << PuzaDEF::Instance()->GE_Polya_Data[k] << " ";
    }
    outFile << endl;
  }
  outFile.close();


  logFile << "[Samples]:" << PuzaDEF::Instance()->GE_Samples << endl;
  logFile << "[Dimensions]:" << PuzaDEF::Instance()->GE_Polya_K << endl;
  logFile << "[Precision]:" << PuzaDEF::Instance()->GE_Polya_Alpha << endl;
  logFile << "[Base Measure]:" << endl;
  for(int k=0; k < PuzaDEF::Instance()->GE_Polya_K; k++){
    logFile << k << ":" << PuzaDEF::Instance()->GE_Polya_Base[k] << " ";
  }
  logFile << endl;
  logFile.close();


  delete[] temp_theta;
  delete[] temp_alpha;
  return 0;
}

#endif
